home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / frasrc19.zip / MPMATH_C.C < prev    next >
Text File  |  1995-03-08  |  17KB  |  649 lines
  1. /* MPMath_c.c (C) 1989, Mark C. Peterson, CompuServe [70441,3353]
  2.      All rights reserved.
  3.  
  4.    Code may be used in any program provided the author is credited
  5.      either during program execution or in the documentation.  Source
  6.      code may be distributed only in combination with public domain or
  7.      shareware source code.  Source code may be modified provided the
  8.      copyright notice and this message is left unchanged and all
  9.      modifications are clearly documented.
  10.  
  11.      I would appreciate a copy of any work which incorporates this code,
  12.      however this is optional.
  13.  
  14.      Mark C. Peterson
  15.      405-C Queen St. Suite #181
  16.      Southington, CT 06489
  17.      (203) 276-9721
  18. */
  19.  
  20. #include <stdlib.h>
  21.  
  22. /* This now in prototyp.h */
  23. /*
  24. #include "mpmath.h"
  25. */
  26.  
  27. #include "prototyp.h"
  28.  
  29. #ifndef XFRACT
  30.  
  31. struct MP *MPsub(struct MP x, struct MP y) {
  32.    y.Exp ^= 0x8000;
  33.    return(MPadd(x, y));
  34. }
  35.  
  36. /* added by TW */
  37. struct MP *MPsub086(struct MP x, struct MP y) {
  38.    y.Exp ^= 0x8000;
  39.    return(MPadd086(x, y));
  40. }
  41.  
  42. /* added by TW */
  43. struct MP *MPsub386(struct MP x, struct MP y) {
  44.    y.Exp ^= 0x8000;
  45.    return(MPadd386(x, y));
  46. }
  47.  
  48. struct MP *MPabs(struct MP x) {
  49.    Ans = x;
  50.    Ans.Exp &= 0x7fff;
  51.    return(&Ans);
  52. }
  53.  
  54. struct MPC MPCsqr(struct MPC x) {
  55.    struct MPC z;
  56.  
  57.         z.x = *pMPsub(*pMPmul(x.x, x.x), *pMPmul(x.y, x.y));
  58.         z.y = *pMPmul(x.x, x.y);
  59.         z.y.Exp++;
  60.    return(z);
  61. }
  62.  
  63. struct MP MPCmod(struct MPC x) {
  64.         return(*pMPadd(*pMPmul(x.x, x.x), *pMPmul(x.y, x.y)));
  65. }
  66.  
  67. struct MPC MPCmul(struct MPC x, struct MPC y) {
  68.    struct MPC z;
  69.  
  70.         z.x = *pMPsub(*pMPmul(x.x, y.x), *pMPmul(x.y, y.y));
  71.         z.y = *pMPadd(*pMPmul(x.x, y.y), *pMPmul(x.y, y.x));
  72.    return(z);
  73. }
  74.  
  75. struct MPC MPCdiv(struct MPC x, struct MPC y) {
  76.    struct MP mod;
  77.  
  78.    mod = MPCmod(y);
  79.         y.y.Exp ^= 0x8000;
  80.         y.x = *pMPdiv(y.x, mod);
  81.         y.y = *pMPdiv(y.y, mod);
  82.    return(MPCmul(x, y));
  83. }
  84.  
  85. struct MPC MPCadd(struct MPC x, struct MPC y) {
  86.    struct MPC z;
  87.  
  88.         z.x = *pMPadd(x.x, y.x);
  89.         z.y = *pMPadd(x.y, y.y);
  90.    return(z);
  91. }
  92.  
  93. struct MPC MPCsub(struct MPC x, struct MPC y) {
  94.    struct MPC z;
  95.  
  96.         z.x = *pMPsub(x.x, y.x);
  97.         z.y = *pMPsub(x.y, y.y);
  98.    return(z);
  99. }
  100.  
  101. struct MPC MPCone = { 0x3fff, 0x80000000l, 0, 0l };
  102.  
  103. struct MPC MPCpow(struct MPC x, int exp) {
  104.    struct MPC z;
  105.    struct MPC zz;
  106.  
  107.    if(exp & 1)
  108.       z = x;
  109.    else
  110.       z = MPCone;
  111.    exp >>= 1;
  112.    while(exp) {
  113.                 zz.x = *pMPsub(*pMPmul(x.x, x.x), *pMPmul(x.y, x.y));
  114.                 zz.y = *pMPmul(x.x, x.y);
  115.                 zz.y.Exp++;
  116.       x = zz;
  117.       if(exp & 1) {
  118.                         zz.x = *pMPsub(*pMPmul(z.x, x.x), *pMPmul(z.y, x.y));
  119.                         zz.y = *pMPadd(*pMPmul(z.x, x.y), *pMPmul(z.y, x.x));
  120.          z = zz;
  121.       }
  122.       exp >>= 1;
  123.    }
  124.    return(z);
  125. }
  126.  
  127. int MPCcmp(struct MPC x, struct MPC y) {
  128.    struct MPC z;
  129.  
  130.         if(pMPcmp(x.x, y.x) || pMPcmp(x.y, y.y)) {
  131.                 z.x = MPCmod(x);
  132.                 z.y = MPCmod(y);
  133.                 return(pMPcmp(z.x, z.y));
  134.    }
  135.    else
  136.       return(0);
  137. }
  138.  
  139. _CMPLX MPC2cmplx(struct MPC x) {
  140.    _CMPLX z;
  141.  
  142.         z.x = *pMP2d(x.x);
  143.         z.y = *pMP2d(x.y);
  144.    return(z);
  145. }
  146.  
  147. struct MPC cmplx2MPC(_CMPLX z) {
  148.    struct MPC x;
  149.  
  150.         x.x = *pd2MP(z.x);
  151.         x.y = *pd2MP(z.y);
  152.    return(x);
  153. }
  154.  
  155. /* function pointer versions added by Tim Wegner 12/07/89 */
  156. /* int        (*ppMPcmp)() = MPcmp086; */
  157. int        (*pMPcmp)(struct MP x, struct MP y) = MPcmp086;
  158. struct MP  *(*pMPmul)(struct MP x, struct MP y)= MPmul086;
  159. struct MP  *(*pMPdiv)(struct MP x, struct MP y)= MPdiv086;
  160. struct MP  *(*pMPadd)(struct MP x, struct MP y)= MPadd086;
  161. struct MP  *(*pMPsub)(struct MP x, struct MP y)= MPsub086;
  162. struct MP  *(*pd2MP)(double x)                 = d2MP086 ;
  163. double *(*pMP2d)(struct MP m)                  = MP2d086 ;
  164. /* struct MP  *(*pfg2MP)(long x, int fg)          = fg2MP086; */
  165.  
  166. void setMPfunctions(void) {
  167.    if(cpu == 386)
  168.    {
  169.       pMPmul = MPmul386;
  170.       pMPdiv = MPdiv386;
  171.       pMPadd = MPadd386;
  172.       pMPsub = MPsub386;
  173.       pMPcmp = MPcmp386;
  174.       pd2MP  = d2MP386 ;
  175.       pMP2d  = MP2d386 ;
  176.       /* pfg2MP = fg2MP386; */
  177.    }
  178.    else
  179.    {
  180.       pMPmul = MPmul086;
  181.       pMPdiv = MPdiv086;
  182.       pMPadd = MPadd086;
  183.       pMPsub = MPsub086;
  184.       pMPcmp = MPcmp086;
  185.       pd2MP  = d2MP086 ;
  186.       pMP2d  = MP2d086 ;
  187.       /* pfg2MP = fg2MP086; */
  188.    }
  189. }
  190.  
  191. #endif /* XFRACT */
  192.  
  193. #ifndef sqr
  194. #define sqr(x) ((x)*(x))
  195. #endif
  196.  
  197. _CMPLX ComplexPower(_CMPLX xx, _CMPLX yy) {
  198.    _CMPLX z, cLog, t;
  199.    double e2x, siny, cosy;
  200.  
  201.    /* fixes power bug - if any complaints, backwards compatibility hook
  202.       goes here TIW 3/95 */
  203.    if(debugflag != 94)
  204.       if(xx.x == 0 && xx.y == 0) {
  205.          z.x = z.y = 0.0;
  206.          return(z);
  207.       }
  208.  
  209.    FPUcplxlog(&xx, &cLog);
  210.    FPUcplxmul(&cLog, &yy, &t);
  211.  
  212.    if(fpu == 387)
  213.       FPUcplxexp387(&t, &z);
  214.    else {
  215.       if(t.x < -690)
  216.          e2x = 0;
  217.       else
  218.          e2x = exp(t.x);
  219.       FPUsincos(&t.y, &siny, &cosy);
  220.       z.x = e2x * cosy;
  221.       z.y = e2x * siny;
  222.    }
  223.    return(z);
  224. }
  225.  
  226. /* 
  227.  
  228.   The following Complex function routines added by Tim Wegner November 1994.
  229.   
  230. */  
  231.  
  232. #define Sqrtz(z,rz) (*(rz) = ComplexSqrtFloat((z).x, (z).y))
  233.  
  234. /* rz=Arcsin(z)=-i*Log{i*z+sqrt(1-z*z)} */
  235. void Arcsinz(_CMPLX z,_CMPLX *rz)
  236. {
  237.   _CMPLX tempz1,tempz2;
  238.  
  239.   if(z.y == 0 && z.x > 1)
  240.   {
  241.      rz->x = HUGE_VAL;
  242.      rz->y = 0;
  243.      return;
  244.   }  
  245.   FPUcplxmul( &z, &z, &tempz1);
  246.   tempz1.x = 1 - tempz1.x; tempz1.y = -tempz1.y;  /* tempz1 = 1 - tempz1 */
  247.   Sqrtz( tempz1, &tempz1);
  248.  
  249.   tempz2.x = -z.y; tempz2.y = z.x;                /* tempz2 = i*z  */  
  250.   tempz1.x += tempz2.x;  tempz1.y += tempz2.y;    /* tempz1 += tempz2 */
  251.   FPUcplxlog( &tempz1, &tempz1);
  252.   rz->x = tempz1.y;  rz->y = -tempz1.x;           /* rz = (-i)*tempz1 */
  253. }   /* end. Arcsinz */
  254.  
  255.  
  256. /* rz=Arccos(z)=-i*Log{z+sqrt(z*z-1)} */  
  257. void Arccosz(_CMPLX z,_CMPLX *rz)
  258. {
  259.   _CMPLX temp;
  260.  
  261.   FPUcplxmul( &z, &z, &temp);
  262.   temp.x -= 1;                                 /* temp = temp - 1 */
  263.   Sqrtz( temp, &temp);
  264.  
  265.   temp.x += z.x; temp.y += z.y;                /* temp = z + temp */
  266.   
  267.   FPUcplxlog( &temp, &temp);
  268.   rz->x = temp.y;  rz->y = -temp.x;              /* rz = (-i)*tempz1 */
  269. }   /* end. Arccosz */
  270.  
  271. void Arcsinhz(_CMPLX z,_CMPLX *rz)
  272. {
  273.   _CMPLX temp;
  274.  
  275.   FPUcplxmul( &z, &z, &temp);
  276.   temp.x += 1;                                 /* temp = temp + 1 */
  277.   Sqrtz( temp, &temp);
  278.   temp.x += z.x; temp.y += z.y;                /* temp = z + temp */
  279.   FPUcplxlog( &temp, rz);
  280. }  /* end. Arcsinhz */
  281.  
  282. /* rz=Arccosh(z)=Log(z+sqrt(z*z-1)} */
  283. void Arccoshz(_CMPLX z,_CMPLX *rz)
  284. {
  285.   _CMPLX tempz;
  286.   FPUcplxmul( &z, &z, &tempz);
  287.   tempz.x -= 1;                              /* tempz = tempz - 1 */
  288.   Sqrtz( tempz, &tempz);
  289.   tempz.x = z.x + tempz.x; tempz.y = z.y + tempz.y;  /* tempz = z + tempz */
  290.   FPUcplxlog( &tempz, rz);
  291. }   /* end. Arccoshz */
  292.  
  293. /* rz=Arctanh(z)=1/2*Log{(1+z)/(1-z)} */
  294. void Arctanhz(_CMPLX z,_CMPLX *rz)
  295. {
  296.   _CMPLX temp0,temp1,temp2;
  297.    
  298.   if( z.x == 0.0){
  299.     rz->x = 0;
  300.     rz->y = atan( z.y);
  301.     return;
  302.   }
  303.   else{
  304.     if( fabs(z.x) == 1.0 && z.y == 0.0){
  305.       return;
  306.     }
  307.     else if( fabs( z.x) < 1.0 && z.y == 0.0){
  308.       rz->x = log((1+z.x)/(1-z.x))/2;
  309.       rz->y = 0;
  310.       return;
  311.     } 
  312.     else{
  313.       temp0.x = 1 + z.x; temp0.y = z.y;             /* temp0 = 1 + z */
  314.       temp1.x = 1 - z.x; temp1.y = -z.y;            /* temp1 = 1 - z */
  315.       FPUcplxdiv( &temp0, &temp1, &temp2);
  316.       FPUcplxlog( &temp2, &temp2);
  317.       rz->x = .5*temp2.x; rz->y = .5*temp2.y;       /* rz = .5*temp2 */
  318.       return;
  319.     }
  320.   }
  321. }   /* end. Arctanhz */
  322.  
  323. /* rz=Arctan(z)=i/2*Log{(1-i*z)/(1+i*z)} */
  324. void Arctanz(_CMPLX z,_CMPLX *rz)
  325. {
  326.   _CMPLX temp0,temp1,temp2,temp3;
  327.   if( z.x == 0.0 && z.y == 0.0)
  328.     rz->x = rz->y = 0;
  329.   else if( z.x != 0.0 && z.y == 0.0){
  330.     rz->x = atan( z.x);
  331.     rz->y = 0;
  332.   }
  333.   else if( z.x == 0.0 && z.y != 0.0){
  334.     temp0.x = z.y;  temp0.y = 0.0;
  335.     Arctanhz( temp0, &temp0);
  336.     rz->x = -temp0.y; rz->y = temp0.x;              /* i*temp0 */
  337.   } 
  338.   else if( z.x != 0.0 && z.y != 0.0){
  339.  
  340.     temp0.x = -z.y; temp0.y = z.x;                  /* i*z */
  341.     temp1.x = 1 - temp0.x; temp1.y = -temp0.y;      /* temp1 = 1 - temp0 */
  342.     temp2.x = 1 + temp0.x; temp2.y = temp0.y;       /* temp2 = 1 + temp0 */
  343.  
  344.     FPUcplxdiv( &temp1, &temp2, &temp3);
  345.     FPUcplxlog( &temp3, &temp3);
  346.     rz->x = -temp3.y*.5; rz->y = .5*temp3.x;           /* .5*i*temp0 */
  347.   }
  348. }   /* end. Arctanz */
  349.  
  350. #define SinCosFudge 0x10000L
  351. #ifdef LONGSQRT
  352. long lsqrt(long f)
  353. {
  354.     int N;
  355.     unsigned long y0, z;
  356.     static long a=0, b=0, c=0;                  /* constant factors */
  357.  
  358.     if (f == 0)
  359.     return f;
  360.     if (f <  0)
  361.     return 0;
  362.  
  363.     if (a==0)                                   /* one-time compute consts */
  364.     {
  365.     a = (long)(fudge * .41731);
  366.     b = (long)(fudge * .59016);
  367.     c = (long)(fudge * .7071067811);
  368.     }
  369.  
  370.     N  = 0;
  371.     while (f & 0xff000000L)                     /* shift arg f into the */
  372.     {                                           /* range: 0.5 <= f < 1  */
  373.     N++;
  374.     f /= 2;
  375.     }
  376.     while (!(f & 0xff800000L))
  377.     {
  378.     N--;
  379.     f *= 2;
  380.     }
  381.  
  382.     y0 = a + multiply(b, f,  bitshift);         /* Newton's approximation */
  383.  
  384.     z  = y0 + divide (f, y0, bitshift);
  385.     y0 = (z>>2) + divide(f, z,  bitshift);
  386.  
  387.     if (N % 2)
  388.     {
  389.     N++;
  390.     y0 = multiply(c,y0, bitshift);
  391.     }
  392.     N /= 2;
  393.     if (N >= 0)
  394.     return y0 <<  N;                        /* correct for shift above */
  395.     else
  396.     return y0 >> -N;
  397. }
  398. #endif
  399. LCMPLX ComplexSqrtLong(long x, long y)
  400. {
  401.    double    mag, theta;
  402.    long      maglong, thetalong;
  403.    LCMPLX    result;
  404.  
  405. #ifndef LONGSQRT
  406.    mag       = sqrt(sqrt(((double) multiply(x,x,bitshift))/fudge +
  407.              ((double) multiply(y,y,bitshift))/ fudge));
  408.    maglong   = (long)(mag * fudge);
  409. #else
  410.    maglong   = lsqrt(lsqrt(multiply(x,x,bitshift)+multiply(y,y,bitshift)));
  411. #endif
  412.    theta     = atan2((double) y/fudge, (double) x/fudge)/2;
  413.    thetalong = (long)(theta * SinCosFudge);
  414.    SinCos086(thetalong, &result.y, &result.x);
  415.    result.x  = multiply(result.x << (bitshift - 16), maglong, bitshift);
  416.    result.y  = multiply(result.y << (bitshift - 16), maglong, bitshift);
  417.    return result;
  418. }
  419.  
  420. _CMPLX ComplexSqrtFloat(double x, double y)
  421. {
  422.    double mag;
  423.    double theta;
  424.    _CMPLX  result;
  425.  
  426.    if(x == 0.0 && y == 0.0)
  427.       result.x = result.y = 0.0;
  428.    else
  429.    {   
  430.       mag   = sqrt(sqrt(x*x + y*y));
  431.       theta = atan2(y, x) / 2;
  432.       FPUsincos(&theta, &result.y, &result.x);
  433.       result.x *= mag;
  434.       result.y *= mag;
  435.    }
  436.    return result;
  437. }
  438.  
  439.  
  440. /***** FRACTINT specific routines and variables *****/
  441.  
  442. #ifndef TESTING_MATH
  443.  
  444. BYTE far *LogTable = (BYTE far *)0;
  445. int MaxLTSize;
  446.  
  447.    /* int LogFlag;
  448.       LogFlag == 1  -- standard log palettes
  449.       LogFlag == -1 -- 'old' log palettes
  450.       LogFlag >  1  -- compress counts < LogFlag into color #1
  451.       LogFlag < -1  -- use quadratic palettes based on square roots && compress
  452.    */
  453.  
  454. void SetupLogTable(void) {
  455.    float l, f, c, m;
  456.    unsigned n, prev, limit, lf;
  457.  
  458.    if (LogFlag > -2) {
  459.       lf = (LogFlag > 1) ? LogFlag : 0;
  460.       if (lf >= (unsigned int)MaxLTSize)
  461.          lf = MaxLTSize - 1;
  462.       Fg2Float((long)(MaxLTSize-lf), 0, m);
  463.       fLog14(m, m);
  464.       Fg2Float((long)(colors-(lf?2:1)), 0, c);
  465.       fDiv(m, c, m);
  466.       for (prev = 1; prev <= lf; prev++)
  467.          LogTable[prev] = 1;
  468.       for (n = (lf?2:1); n < (unsigned int)colors; n++) {
  469.          Fg2Float((long)n, 0, f);
  470.          fMul16(f, m, f);
  471.          fExp14(f, l);
  472.          limit = (unsigned int)Float2Fg(l, 0) + lf;
  473.          if (limit > (unsigned int)MaxLTSize || n == (unsigned int)(colors-1))
  474.             limit = MaxLTSize;
  475.          while (prev <= limit)
  476.             LogTable[prev++] = (BYTE)n;
  477.       }
  478.    } else {
  479.       if ((lf = 0 - LogFlag) >= (unsigned int)MaxLTSize)
  480.          lf = MaxLTSize - 1;
  481.       Fg2Float((long)(MaxLTSize-lf), 0, m);
  482.       fSqrt14(m, m);
  483.       Fg2Float((long)(colors-2), 0, c);
  484.       fDiv(m, c, m);
  485.       for (prev = 1; prev <= lf; prev++)
  486.          LogTable[prev] = 1;
  487.       for (n = 2; n < (unsigned int)colors; n++) {
  488.          Fg2Float((long)n, 0, f);
  489.          fMul16(f, m, f);
  490.          fMul16(f, f, l);
  491.          limit = (unsigned int)(Float2Fg(l, 0) + lf);
  492.          if (limit > (unsigned int)MaxLTSize || n == (unsigned int)(colors-1))
  493.             limit = MaxLTSize;
  494.          while (prev <= limit)
  495.             LogTable[prev++] = (BYTE)n;
  496.       }
  497.    }
  498.    LogTable[0] = 0;
  499.    if (LogFlag != -1)
  500.       for (n = 1; n < (unsigned int)MaxLTSize; n++) /* spread top to incl unused colors */
  501.          if (LogTable[n] > LogTable[n-1])
  502.             LogTable[n] = (BYTE)(LogTable[n-1]+1);
  503. }
  504.  
  505. long far ExpFloat14(long xx) {
  506.    static float fLogTwo = (float)0.6931472;
  507.    int f;
  508.    long Ans;
  509.  
  510.    f = 23 - (int)RegFloat2Fg(RegDivFloat(xx, *(long*)&fLogTwo), 0);
  511.    Ans = ExpFudged(RegFloat2Fg(xx, 16), f);
  512.    return(RegFg2Float(Ans, (char)f));
  513. }
  514.  
  515. double TwoPi;
  516. _CMPLX temp, BaseLog;
  517. _CMPLX cdegree = { 3.0, 0.0 }, croot   = { 1.0, 0.0 };
  518.  
  519. int ComplexNewtonSetup(void) {
  520.    threshold = .001;
  521.    periodicitycheck = 0;
  522.    if(param[0] != 0.0 || param[1] != 0.0 || param[2] != 0.0 ||
  523.       param[3] != 0.0) {
  524.       croot.x = param[2];
  525.       croot.y = param[3];
  526.       cdegree.x = param[0];
  527.       cdegree.y = param[1];
  528.       FPUcplxlog(&croot, &BaseLog);
  529.       TwoPi = asin(1.0) * 4;
  530.    }
  531.    return(1);
  532. }
  533.  
  534. int ComplexNewton(void) {
  535.    _CMPLX cd1;
  536.  
  537.    /* new = ((cdegree-1) * old**cdegree) + croot
  538.             ----------------------------------
  539.                  cdegree * old**(cdegree-1)         */
  540.  
  541.    cd1.x = cdegree.x - 1.0;
  542.    cd1.y = cdegree.y;
  543.  
  544.    temp = ComplexPower(old, cd1);
  545.    FPUcplxmul(&temp, &old, &new);
  546.  
  547.    tmp.x = new.x - croot.x;
  548.    tmp.y = new.y - croot.y;
  549.    if((sqr(tmp.x) + sqr(tmp.y)) < threshold)
  550.       return(1);
  551.  
  552.    FPUcplxmul(&new, &cd1, &tmp);
  553.    tmp.x += croot.x;
  554.    tmp.y += croot.y;
  555.  
  556.    FPUcplxmul(&temp, &cdegree, &cd1);
  557.    FPUcplxdiv(&tmp, &cd1, &old);
  558.    if(DivideOverflow)
  559.    {
  560.       DivideOverflow = 0;
  561.       return(1);
  562.    }
  563.    new = old;
  564.    return(0);
  565. }
  566.  
  567. int ComplexBasin(void) {
  568.    _CMPLX cd1;
  569.    double mod;
  570.  
  571.    /* new = ((cdegree-1) * old**cdegree) + croot
  572.             ----------------------------------
  573.                  cdegree * old**(cdegree-1)         */
  574.  
  575.    cd1.x = cdegree.x - 1.0;
  576.    cd1.y = cdegree.y;
  577.  
  578.    temp = ComplexPower(old, cd1);
  579.    FPUcplxmul(&temp, &old, &new);
  580.  
  581.    tmp.x = new.x - croot.x;
  582.    tmp.y = new.y - croot.y;
  583.    if((sqr(tmp.x) + sqr(tmp.y)) < threshold) {
  584.       if(fabs(old.y) < .01)
  585.          old.y = 0.0;
  586.       FPUcplxlog(&old, &temp);
  587.       FPUcplxmul(&temp, &cdegree, &tmp);
  588.       mod = tmp.y/TwoPi;
  589.       coloriter = (long)mod;
  590.       if(fabs(mod - coloriter) > 0.5) {
  591.          if(mod < 0.0)
  592.             coloriter--;
  593.          else
  594.             coloriter++;
  595.       }
  596.       coloriter += 2;
  597.       if(coloriter < 0)
  598.          coloriter += 128;
  599.       return(1);
  600.    }
  601.  
  602.    FPUcplxmul(&new, &cd1, &tmp);
  603.    tmp.x += croot.x;
  604.    tmp.y += croot.y;
  605.  
  606.    FPUcplxmul(&temp, &cdegree, &cd1);
  607.    FPUcplxdiv(&tmp, &cd1, &old);
  608.    if(DivideOverflow)
  609.    {
  610.       DivideOverflow = 0;
  611.       return(1);
  612.    }
  613.    new = old;
  614.    return(0);
  615. }
  616.  
  617. /*
  618.  * Generate a gaussian distributed number.
  619.  * The right half of the distribution is folded onto the lower half.
  620.  * That is, the curve slopes up to the peak and then drops to 0.
  621.  * The larger slope is, the smaller the standard deviation.
  622.  * The values vary from 0+offset to range+offset, with the peak
  623.  * at range+offset.
  624.  * To make this more complicated, you only have a
  625.  * 1 in Distribution*(1-Probability/Range*con)+1 chance of getting a
  626.  * Gaussian; otherwise you just get offset.
  627.  */
  628. int GausianNumber(int Probability, int Range) {
  629.    int n, r;
  630.    long Accum = 0, p;
  631.  
  632.    p = divide((long)Probability << 16, (long)Range << 16, 16);
  633.    p = multiply(p, con, 16);
  634.    p = multiply((long)Distribution << 16, p, 16);
  635.    if(!(rand15() % (Distribution - (int)(p >> 16) + 1))) {
  636.       for(n = 0; n < Slope; n++)
  637.          Accum += rand15();
  638.       Accum /= Slope;
  639.       r = (int)(multiply((long)Range << 15, Accum, 15) >> 14);
  640.       r = r - Range;
  641.       if(r < 0)
  642.          r = -r;
  643.       return(Range - r + Offset);
  644.    }
  645.    return(Offset);
  646. }
  647.  
  648. #endif
  649.